#include "header.h"

/* ----------------------------------------------------------------------- */
main(int argc, char* argv[]){
/* ----------------------------------------------------------------------- */
int i;
struct node *p;
FILE* ferr;
void* nullpt;

void title();
void read_grid(struct node**);
void connstat(struct node*);
void svd_test(struct node*);
void poly_test(struct node*);
void initialize_var(struct node*);
void ShowConn(struct node*);
void start(struct node*);


if(argc == 1){ printf("Specify input file\n"); exit(0); }
InputFile = argv[1];
title();
read_input(); printf("Reading input data from %s\n",InputFile);
read_grid(&p);
connstat(p);
geometric(p);
svd_test(p);
poly_test(p);
initialize_var(p);
parameters_used();
printf("Please see the file param.dat for the input values.\n");
//ShowConn(p); exit(0);
start(p);
output(p);
printf("Results written into output.dat\n");
}

/* ----------------------------------------------------------------------- */
void read_input(){
/* ----------------------------------------------------------------------- */
int i;
char dummy[50];
FILE* fptr;

fptr = fopen(InputFile,"r");
if(fptr == NULL){
	printf("Input file is missing.\n");
	exit(0);
	}
fscanf(fptr,"%s %lf",dummy,&mach_inf);
fscanf(fptr,"%s %lf",dummy,&alpha_inf);
fscanf(fptr,"%s %lf",dummy,&beta_inf);
fscanf(fptr,"%s %lf",dummy,&cfl);
fscanf(fptr,"%s %lf",dummy,&wall_cfl);
fscanf(fptr,"%s %lf",dummy,&swbc_relax);
fscanf(fptr,"%s %lf",dummy,&KKK);
fscanf(fptr,"%s %d",dummy,&rk_stage);
fscanf(fptr,"%s %d",dummy,&NOIMP);
fscanf(fptr,"%s %lf",dummy,&MIN_RESIDUE);
fscanf(fptr,"%s %d",dummy,&LAST_ITER);
fscanf(fptr,"%s %d",dummy,&MAX_ITER);
fscanf(fptr,"%s %d",dummy,&ResFreq);
fscanf(fptr,"%s %d",dummy,&ResDisp);
fscanf(fptr,"%s %d",dummy,&startmode);
fscanf(fptr,"%s %d",dummy,&time_mode);
fscanf(fptr,"%s %d",dummy,&qmode);
fscanf(fptr,"%s %d",dummy,&wallmode);
fscanf(fptr,"%s %d",dummy,&no_spl_points);
for(i=0;i<no_spl_points;i++) fscanf(fptr,"%d",&spl_point[i]);
fscanf(fptr,"%s",GridFile);
fscanf(fptr,"%s",ParFile);
fscanf(fptr,"%s",ResFile);
fscanf(fptr,"%s",OutFile);
fscanf(fptr,"%s",PointFile);
fscanf(fptr,"%s",VigieFile);
fscanf(fptr,"%s",CpFile);
fclose(fptr);

alpha_inf = M_PI*alpha_inf/180.0;
beta_inf  = M_PI*beta_inf/180.0;
}
/* ------------------------------------------------------------------------ */
void parameters_used(){
/* ------------------------------------------------------------------------ */
FILE* fptr;

fptr = fopen(ParFile,"w");
fprintf(fptr,"-----------------------------------------\n");
fprintf(fptr,"    Input parameters for 2D KMM Code     \n");
fprintf(fptr,"-----------------------------------------\n");
fprintf(fptr,"No. of grid points = %d\n",no_grid_points);
fprintf(fptr,"Min edge length    = %e\n",MinEdge);
fprintf(fptr,"Max edge length    = %e\n",MaxEdge);
fprintf(fptr,"Length for venkat  = %e\n",DsVenkat);
fprintf(fptr,"Max condition no.  = %f, %f\n",MaxCond,MaxCondWeighted);
fprintf(fptr,"Mach number        = %f\n",mach_inf);
fprintf(fptr,"Angle of attack    = %f\n",alpha_inf*180.0/M_PI);
fprintf(fptr,"Yaw angle          = %f\n",beta_inf*180.0/M_PI);
fprintf(fptr,"Temperature        = %f\n",temp_inf);
fprintf(fptr,"Density            = %f\n",rho_inf);
fprintf(fptr,"CFL number         = %f\n",cfl);
fprintf(fptr,"Wall CFL number    = %f\n",wall_cfl);
fprintf(fptr,"SWBC relaxation    = %f\n",swbc_relax);
fprintf(fptr,"Venkat             = %f\n",KKK);
fprintf(fptr,"No. of RK stages   = %d\n",rk_stage);
fprintf(fptr,"No. of impli iter  = %d\n",NOIMP);
fprintf(fptr,"Min. Residue       = %.2e\n",MIN_RESIDUE);
fprintf(fptr,"Max. no. of iter   = %d\n",MAX_ITER);

if(qmode == QLIMIT) 
	fprintf(fptr,"q Derivative       = Limited\n");
else 
	fprintf(fptr,"q Derivative       = Unlimted\n");

if(wallmode == NWBC)
	fprintf(fptr,"Wall treatement    = NWBC with correction\n");
else
	fprintf(fptr,"Wall treatement    = NWBC\n");

if(time_mode == GLOBAL)
	fprintf(fptr,"Time - stepping    = Global\n");
else
	fprintf(fptr,"Time - stepping    = Local\n");
fprintf(fptr,"Grid file          = %s\n",GridFile);
fclose(fptr);
}
/* ------------------------------------------------------------------------ */
void initialize_var(struct node* p){
/* ------------------------------------------------------------------------ */
int i;
double velocity,sonic,dist,rho,ru,rv,rw,E,fdum;
FILE* fpt;

rho_inf  = 1.0;
velocity = 1.0;

u_inf 	 = velocity*cos(alpha_inf);
w_inf 	 = 0.0;
v_inf 	 = velocity*sin(alpha_inf);

/*
u_inf 	 = velocity*cos(alpha_inf);
v_inf 	 = 0.0;
w_inf 	 = velocity*sin(alpha_inf);
*/

sonic	   = velocity/mach_inf;
temp_inf = sonic*sonic/(GAMMA*GAS_CONST);
q_inf    = 0.5*rho_inf*velocity;

pre_inf  = rho_inf*GAS_CONST*temp_inf;
H_inf    = sonic*sonic/(GAMMA-1.0) + 0.5*velocity*velocity;
S_inf    = pre_inf/pow(rho_inf, GAMMA);


switch(startmode){
	case SCRATCH:
		printf("Initializing all variables to free-stream values ....... \n");
		for(i=0;i<no_grid_points;i++){
			p[i].var[0] = rho_inf;
			p[i].var[1] = u_inf;
			p[i].var[2] = v_inf;
			p[i].var[3] = w_inf;
			p[i].var[4] = pre_inf;

			}
		break;
	case RESTART:
		printf("Starting from previous output file ....... \n");
		fpt = fopen("INIT","r");
		for(i=0;i<no_grid_points;i++){
			fscanf(fpt,"%lf%lf%lf%lf%lf%lf%lf%lf",&fdum,&fdum,&fdum,&rho,
                &ru,&rv,&rw,&E);
			if(rho<=0.0){
				printf("Restart: Negative density.\n");
				exit(0);
				}
			if(E<=0.0){
				printf("Restart: Negative Energy density.\n");
				exit(0);
				}
			p[i].var[0] = rho;
			p[i].var[1] = ru;
			p[i].var[2] = rv;
			p[i].var[3] = rw;
			p[i].var[4] = E;
			if(p[i].var[4]<=0.0){
				printf("Restart: Negative pressure .\n");
				printf("Node = %d, x = %f, y = %f, z = %f\n",
                   i, p[i].x, p[i].y, p[i].z);
				exit(0);
				}
			}
		break;
	}
printf("Free stream velocity = %f, %f, %f\n", u_inf, v_inf, w_inf); 
}

/* ------------------------------------------------------------------------ */
